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Frustrated magnetism on the kagome lattice has been a fertile ground for rich and fascinating 
physics, ranging from experimental evidence of a spin liquid to theoretical predictions of exotic 
superconductivity. Among experimentally realized spin-| kagome magnets, herbertsmithite, kapel- 
lasite, and haydeeite [(Zn,Mg)Cu 3 (OH) 6 Cl 2 ] are all well described by a three-parameter Heisenberg 
model, but they exhibit distinctly different physics. We address the problem using a pseudofermion 
functional renormalization-group approach and analyze the low-energy physics in the experimen¬ 
tally accessible parameter range. Our analysis places kapellasite and haydeeite near the boundaries 
between magnetically ordered and disordered phases, implying that slight modifications could dra¬ 
matically affect their magnetic properties. Inspired by this, we perform ab initio density functional 
theory calculations of (Zn,Mg,Cd)Cu 3 (OH) 6 Cl 2 at various pressures. Our results suggest that by 
varying pressure and composition one can traverse a paramagnetic regime between different mag¬ 
netically ordered phases. 

PACS numbers: 75.10.Jm, 75.10.Kt, 75.40.Mg, 05.10.Cc 


Introduction. Quantum magnetism in low-dimensional 
systems with parametric or geometrical frustration has 
been a highly inspiring field of research ever since the 
seminal paper of Pomeranchuk [1] . A Holy Grail of the 
field has been the spin-^ antiferromagnet on the kagome 
lattice (Fig. 1), where the geometrical frustration inher¬ 
ent in the individual triangles is only marginally allevi¬ 
ated through a corner sharing lattice pattern [2]. It is 
widely conjectured to host a spin liquid phase, albeit the 
nature and the topological classification of this phase are 
still controversial [3-15]. It may host exotic supercon¬ 
ducting phases [16, 17]. 

Recently, significant progress has been achieved re¬ 
garding experimental realizations of kagome magnets. 
The couplings in the actual materials, however, differ 
significantly from idealized models. The most prominent 
materials are herbertsmithite [ZnGu 3 (OH) 6 Gl 2 ] [18, 19] 
and its polymorphs, kapellasite [ZnCus (OH)gCl 2 ] [20, 
21] and haydeeite [MgGu 3 (OH) 6 Gl 2 ] [22-26], with 
ground states ranging from potentially paramagnetic 
(PM) to weakly ferromagnetic (FM) phases. This variety 
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of ground states in structurally similar systems calls for a 
thorough theoretical investigation of the experimentally 
relevant couplings and their implications. In this Rapid 
Communication, we focus on the following aspects: (1) 
Which phases—ordered or PM—are realized, depending 
on Heisenberg couplings Ji, J 2 , and Jd (Fig. 1)? (2) 
What determines the nature of magnetic interactions, 
and why are they so different in these compounds? (3) 
Can we deliberately manipulate these materials to probe 
different parts of the phase diagram, and, in particular, 
the PM (possibly spin-liquid) phase indicated in Fig. 2? 

To begin with, we map out the zero-temperature phase 
diagram of the Ji-J 2 -Jd kagome model. Specifically, 
we employ the pseudofermion functional renormalization 
group (PF-FRG) [27-29] method to compute magnetic 
fluctuations. We find that the experimentally estimated 
parameters place kapellasite near the borderline between 
the PM and the antiferromagnetically ordered cuboc-2 
phase and haydeeite near the PM/FM border. Using 
ab initio density functional theory (DFT) calculations, 
we then discuss the reliability of these parameters and 
possible microscopic origins for their variations. We pro¬ 
ceed with suggestions on how one can modify these com¬ 
pounds to shift them away from their current positions 
and explore other parts of the phase diagram. Among 
other aspects, our results provide an independent assess¬ 
ment of the initial placement of the materials in the phase 
diagram. 

Herbertsmithite, kapellasite (KL), and haydeeite (HD) 
feature geometrically perfect Cu^+ -S' = ^ kagome planes 
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FIG. 1. For all kagome magnets considered, the lattice is 
formed by 00 ^+ S=\ spins (black), and the Heisenberg ex¬ 
change couplings are given by nearest neighbor Ji, second 
nearest neighbor J 2 , and diagonal Jd across the hexagons. 
Kapellasite and haydeeite feature in-plane Zn^+ and Mg^+ 
ions (green), respectively, at the center of the hexagons, trig¬ 
gering appreciable values for Jd- 


with the nearest-neighbor (NN) superexchange Ji medi¬ 
ated by 0H“ and Cl“. The minimal model also includes 
subdominant interactions J 2 and J 3 and significant Jd 
(Fig. 1). Experiment [30] and calculations [31, 32] sug¬ 
gest that in herbertsmithite, where only Cu is present in 
the kagome planes, both J 3 and Jd are negligible, and 
thus the material is well described by a NN antiferro¬ 
magnetic model with Ji « 180 K, with a small but non- 
negligible J 2 [29, 33]. The quantum paramagnetic ground 
state of such an antiferromagnet has been intensively 
studied theoretically (see, e.g.. Refs. [3-9, 11-15, 34]), 
and there is experimental evidence of a spin liquid state 
in herbertsmithite [19, 35, 36]. 

On the other hand, in KL and HD, the Zn and Mg ions, 
respectively, occupy the centers of the Cu hexagons [21], 
thus spanning the Cu pairs connected by Jd- This seems 
to explain why the Jd interaction is sizable, and differs 
considerably between the two compounds. However, the 
nature of this interaction is probably more complex than 
that, as discussed in detail in the Supplemental Mate¬ 
rial [37]. Bernu et aL [38] extracted J’s in KL from 
the temperature dependencies of magnetic susceptibil¬ 
ity and specific heat, while Boldrin et aL [26] did 
the same for HD, using spin-wave dispersion. With the 
caveat that these are distinctly different experimental 
procedures, the estimated exchange coupling constants 
are (Ji,J 2 ,Jd) = (-12,-4,15.6) K for KL [38] and 
(Ji, Jrf) = (-38,11) K, with J 2 /J 1 < 0.1 for HD [26]. 
The small and negative values of Ji signal a large cancel¬ 
lation of the anti- and ferromagnetic contributions to the 
NN superexchange, which, as we discuss in the Supple¬ 
mental Material [37], is quite unexpected, but that they 
are close in both compounds is consistent with their simi¬ 
lar geometries. Further investigations using electron spin 


resonance estimated the symmetric exchange anisotropy 
D to be only |D/Ji| ^ 3% [39] , thus justifying the use 
of the {Ji,J 2 ,Jd) isotropic Heisenberg model as a good 
starting point for both KL and HD. 

Model and Methods- The Heisenberg Hamiltonian 
reads 

H = Ji Vs,-S, + J2 Vs,-s, + Jd V s,-s„ (1) 

(b> ((b)> (((b»>d 

where is the spin-^ operator at site i- Here, Ji, J 2 < 0 
(ferromagnetic) and Jdi^ 0 (antiferromagnetic), normal¬ 
ized so that I Ji| -I- IJ 2 I + Jd = I- The indices {ij), {{ij)), 
and {{{ij)))d denote sums over NN and next-NN bonds, 
and the diagonals of hexagons, respectively (Fig. 1). 

In the PF-FRG approach [27-29, 40, 41], we first 
rewrite Eq. (1) in terms of pseudofermions as = 

I J2ap {a, P =t, 4 .), where Ci^a are the pseud¬ 

ofermion operators, and cr are Pauli matrices. This en¬ 
ables us to apply Wick’s theorem and develop a diagram¬ 
matic technique. To this end, an infrared frequency cut¬ 
off A is introduced in the fermion propagator. The FRG 
ansatz (for recent reviews, see, e.g.. Refs. [42, 43]) then 
formulates an infinite hierarchy of coupled integrodiffer- 
ential equations for the evolution of all m-particle vertex 
functions under the flow of A. Within PF-FRG, the trun¬ 
cation of this system of equations to a closed set is ac¬ 
complished by the inclusion of only two-particle reducible 
two-loop contributions, which ensures sufficient backfeed- 
ing of the self-energy corrections to the two-particle ver¬ 
tex evolution [44]. A crucial advantage of the PF-FRG 
is that the diagrammatic summation incorporates vertex 
corrections between all interaction channels, i.e., treats 
magnetic ordering and disordering tendencies on equal 
footing. The PF-FRG equations are solved numerically 
by discretizing the frequency dependencies of the vertex 
functions and limiting the spatial dependencies to a finite 
cluster. We used 64 discretized frequencies and a clus¬ 
ter of 432 sites. In the PF-FRG approach, the onset of 
magnetic long-range order is signaled by a breakdown of 
the smooth RG flow, whereas a smooth evolution down 
to A —> 0 (where A is the infrared frequency cutoff) in¬ 
dicates PM behavior [27] (Fig. 3). From the effective 
low-energy two-particle vertex, we obtain the spin-spin 
correlation function in real space, which we then con¬ 
vert into the momentum-resolved spin susceptibility (see 
Fig. 4). 

Previous applications of the PF-FRG method to frus¬ 
trated magnets have been extremely successful. In par¬ 
ticular, (i) the determined magnetic and nonmagnetic 
phases of the spin-^ Heisenberg J1-J2 antiferromagnet 
on the square lattice and the locations of the phase tran¬ 
sitions quantitatively agree with DMRG, exact diagonal- 
ization, and other methods [27], (ii) the phase diagram 
of the J1-J2-J3 Heisenberg model on the honeycomb lat¬ 
tice agrees perfectly with exact diagonalization [40, 45], 
(iii) the phase diagram of the Kitaev-Heisenberg model 
is also correctly determined within the PF-FRG [41, 46], 
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FIG. 2. Quantum phase diagram of the Ji-J 2 -Jd Heisenberg model as defined in Eq. (1). It features a large paramagnetic (PM) 
domain for intermediate Jd- The exchange couplings estimated from fitting experimental data for kapellasite and haydeeite 
are marked [26, 38]. The static spin structure factors in the extended Brillouin zone are shown next to each phase. The 
corresponding classical phase diagram is shown in the upper left. The evolution of the couplings in the different materials 
under application of pressure, as calculated by ab initio methods, is shown in the enlarged region on the right. 


in particular, the short range nature of the spin corre¬ 
lations in the Kitaev limit is correctly reproduced, and 
(iv) the spin structure factor of the NN Heisenberg anti- 
ferromagnet on the Kagome lattice in PF-FRG is in very 
good quantitative agreement with DMRG [9, 29], which 
is of particular relevance to the problem at hand. 

For ab initio DFT calculations we used the generalized 
gradient approximation (GGA) functional [47]. Struc¬ 
ture optimizations were performed with the projector 
augmented wave method within the VASP code [48, 49], 
and accurate total energies were calculated using the all¬ 
electron FPLO code [50]. 

Results. The PF-FRG quantum phase diagram of the 
Ji-J 2 -Jd model of Eq. (1) is depicted in Fig. 2. Indi¬ 
vidual data points are labeled according to which type of 
phase they belong to in the PF-FRG. For small Jd, FM 
dominates. For intermediate Jd and large J 2 , a -\/3 x 
order is found which changes into the cuboc-1 order for 
increasing Jd- The cuboc orders describe 12-sublattice 
noncoplanar orders in which the spins orient towards the 
corners of a cuboctahedron [51, 52] (that is, along the 12 
possible [110] directions). For the domains discussed so 
far, the quantum phase diagram approximately matches 
the classical phase diagram of Eq. (1) (Fig. 2). Quantum 
corrections start to become visible closer to the cuboc- 
l/cuboc-2 boundary, as J 2 becomes smaller than Ji (for 


large Jd). The classical first order transition line be¬ 
tween the cuboc-1 and cuboc-2 phases is then replaced 
by a narrow vertical strip (Ji « J 2 ) in the quantum 
case, depicted by a merging gradient in Fig. 2, where an 
effectively ID paramagnetic chain regime is found. As 
the most important modification to the classical picture, 
however, an extended PM regime emerges for J 2 /J 1 < 1 
separating the cuboc-2 from the FM domain. Its spin 
susceptibility profile has a well-defined wave-vector de¬ 
pendence featuring dominant short-distance correlations 
with soft maxima at cuboc-2 ordering wave vectors and 
subdominant FM correlations. This type of magnetic 
fluctuation profile is rather peculiar for a PM phase and 
fundamentally different from what is found for herbert- 
smithite [29, 30]. As we enter the PM regime from the 
cuboc-2 phase by lowering Jd, the magnetic correlations 
change quantitatively but not qualitatively, as is mani¬ 
fest from a comparison of their spin susceptibility pro¬ 
files in Fig. 4(c) and Fig. 4(d). The only notable differ¬ 
ence is more spectral weight smearing in the PM regime. 
Within the PM regime, the spin correlations are found to 
be short-ranged, and calculations of the dimer response 
function rule out any kind of valence-bond crystal or¬ 
der up to a 36-site cell. Note that the Ji-J 2 -Jd model 
has been recently analyzed by variational Monte Carlo 
(VMC) [53] . There, in the large Jd regime, noncoplanar 
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FIG. 3. Representative RG flows of the magnetic susceptibil¬ 
ities at the ordering wave vectors of the four ordered regimes 
of Fig. 2 and the paramagnetic regime. The points at which 
the solid lines become dashed (marked by arrows) indicate an 
instability in the flow and express the onset of order. In the 
smooth flow (red curve) indicating paramagnetism, no such 
instability is found. The small oscillations below A ~ 0.1 in 
this flow are due to frequency discretization. 

cuboc-1 order is absent, and cuboc-2 order is reduced to 
a small part of the parameter space. Instead, for a signif¬ 
icant range of Jd and depending on J 2 /^i, two distinct 
gapless 17(1) chiral spin liquids with a spinon Fermi sur¬ 
face are found over an extended region. As opposed to 
the PF-FRG, which treats magnetic order and disorder 
tendencies on the same footing, a certain bias of VMC 
against cuboc-1 and cuboc-2 orders can be argued for 
on the basis of the variational wave functions employed. 
The noncoplanar structure of cuboc-1 and cuboc-2 orders 
implies that the corresponding chosen Jastrow wave func¬ 
tions are inaccurate as the Jastrow factor does not cor¬ 
rectly describe the relevant quantum fluctuations on top 
of the classical state [53]. On the other hand, PF-FRG 
does not suffer from this deficiency, and, if anything, may 
slightly overestimate the PM domain. 

The location of the high-temperature series expansion 
estimate of exchange couplings for KL [38] is marked 
within the center triangle by a star at (0.38,0.13,0.49) 
in Fig. 2. This is very close to the boundary between 
the cuboc-2 and the quantum PM phases. Experimen¬ 
tally, KL shows no spin freezing and persistent fluctua¬ 
tions down to 20 mK (by /rSR), a diffused continuum of 
excitations (inelastic neutron scattering), and the diver¬ 
gence of the intrinsic local susceptibility for T —)■ 0 in 
NMR [54] . The static spin structure factor shows a well- 
defined wave-vector dependence exhibiting AFM short- 
range correlations [51, 55], consistent with the cuboc-2 
pattern. This whole set of experiments has been inter¬ 
preted in favor of a gapless quantum spin liquid [54] 
close to the cuboc-2 AFM order. The delicate location 
of KL might have important experimental implications, 
in that only moderate modifications in material synthesis 
or experimental conditions amounting to strain, pressure. 


defects/impurities, and the imminent presence of differ¬ 
ent type of anisotropies would lead to significant effects. 
Slight modifications of the Heisenberg coupling constants 
could drive kapellasite either into a weak cuboc-2 order 
or towards a quantum PM phase. As tentatively ob¬ 
served in current experiments, this finding is consistent 
with the compound exhibiting strong magnetic frustra¬ 
tion and significant ordering fluctuation tendencies to¬ 
wards cuboc-2 at the same time. Our PF-FRG calcula¬ 
tions show that this kapellasite location yields a critical 
flow, that is, neither shows a robust and smooth RG flow 
down to A = 0 that would point at quantum PM nor ex¬ 
hibits a clear signature of an order-induced breakdown. 

The location of the linear spin-wave estimate of ex¬ 
change couplings for haydeeite [26] is marked by the di¬ 
amond at (0.77,0.0,0.23) in Fig. 2. Remarkably, it is 
likewise located at the border with the PM regime, but 
now on the FM side. Experiments [26] suggest a very 
weak EM order below 4.2 K. From our PF-FRG analy¬ 
sis, ferromagnetic order is weak but unambiguous, and 
the magnetic fluctuations clearly show FM signatures: a 
dominant peak at F and subdominant peaks at Mg [see 
Fig. 4(e)]. 

To summarize, the reported exchange parameters [26, 
38] place KL and HD on opposite sides of the (arguably 
most interesting) paramagnetic region. Moreover, it 
seems like these two compounds, accidentally, are both at 
or very close to a borderline between two phases. While 
making them especially intriguing, it also considerably 
complicates their study. It is therefore highly desirable 
to be able to modify the same compounds continuously, 
in order to “traverse” the phase diagram. In principle, 
there are several ways to do so. 

One option is alloying Mg and Zn by creating a 
mixed compound Mga;Zni_a;Gu 3 (OH)eGl 2 as suggested 
in Ref. [26] . However, while this proposal creates a sys¬ 
tem with average exchange couplings intermediate be¬ 
tween those in KL and HD, in reality it will have random 
bonds with exchange constants similar to those either in 
KL or HD and is more likely to freeze into a spin glass 
state rather than to develop a spin-liquid phase. Besides, 
chemical substitution of Zn by Mg naturally affects J^, 
but the effect on Ji is harder to predict. 

These considerations lead us to propose alternative op¬ 
tions without introducing additional disorder. For that, 
we need first to exercise some caution when using the 
experimental numbers. Indeed, Refs. [26, 38] are com¬ 
plementary in terms of methodology used to extract the 
exchange coupling constants; Ref. [38] relies on mag- 
netometry/calorimetry while Ref. [26] does a spin-wave 
analysis. In both cases, and especially in HD, actual 
samples have considerable excess of Zn or Mg, substitut¬ 
ing for Cu. Each missing Gu creates four incompletely 
frustrated spins, which may alter the results compared 
to the stoichiometric compound. While the reported pa¬ 
rameters for HD are consistent with a FM ground state, 
the experimental data [26] looks more complicated than 
that. Indeed, the ordered moment from neutron scat- 
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FIG. 4. (a) The first (solid fine) and extended (dashed fine) Brillouin zones of the kagome lattice, (b)-(e) Representative 

spin susceptibility profiles obtained in PF-FRG for different regimes of the quantum phase diagram in Fig. 2: (b) cuboc-1, (c) 
cuboc-2, (d) paramagnetic (PM), (e) ferromagnetic (FM). 


tering was estimated to be < 0.2 ^b [26], while the 

saturation moment from magnetization was 0.83 [26] 
or I.O/tb [25], and the spontaneous ordered moment is 
0.02^b [24] or less [25]. The Curie-Weiss effective mo¬ 
ment is 1.83 ^b [24], consistent with an ordered moment 
greater than I/xb. A 40 times difference between the 
saturation and spontaneous magnetization is highly un¬ 
usual, and so is the discrepancy between spontaneous 
and neutron-measured moments. While this has been 
vaguely ascribed to frustration [26], the “frustration pa¬ 
rameter”, usually defined as the ratio of the mean field 
Tew and frustration-suppressed T^ is in fact less than 
one here. The magnetic susceptibility starts growing with 
cooling below 5 K, but does not diverge at the putative 
Tc = 4.2 K, and instead starts flattening out below 4 K. 
Overall, the magnetic behavior described in Ref. [26] is 
more typical for canted antiferromagnets than for ferro- 
magnets. On the other hand, an independent study [25] 
found much larger magnetic moment and stronger ferro¬ 
magnetic behavior. 

Given the experimental situation, we decided to ad¬ 
dress the question theoretically by performing DFT-|-t7 
calculations (see Supplemental Material). We used the 
FPLO program with [/ = 6 — 8 eV, using the fully- 
localized double counting scheme, orthogonal projection 
of 3d density, and gradient corrections in the DFT func¬ 
tional, and we found that, in agreement with previous 
observations [32], this setup gives the closest agreement 
with the experiment for KL. As discussed in detail in the 
Supplemental Material, the exchange coupling constants 
in these materials are not only small, but they also de¬ 
pend on the technical details of the setup (we remind 
the reader that while DFT is a first principles method, 
DFT-|-[/ is not), which accounts, for instance, for the dif¬ 
ference between Refs. [32] and [55]. On the other hand, 
the trends in the dependence of the exchange parameters 
on the geometry and chemistry in these compounds are 
quite robust, and therefore such calculations can be used 
as a guidance for modifying existing materials in order 
to steer them toward one or another magnetic phase. 

With this setup, using reported crystal structures (and 
optimized positions of hydrogen, since these are not 
known experimentally) we obtained for KL (Ji, J 2 , J^) = 
(-12.5,-0.55,16.1)K, corresponding to (0.43,0.02,0.55) 


. This is rather close to (—12, —4,15.6)K [38], albeit 
a bit deeper in the cuboc-2 phase. For HD we find 
(—21.2,0.57,12.8)K, or (0.61,0.02, 0.37), placing it in the 
PM regime (see inset in Fig. 2). By comparing calcula¬ 
tions for the same structure, but substituting Zn for Mg, 
or for the same composition but different structure, we 
found, not surprisingly, that Ji is predominantly (80- 
90%) controlled by the structure and Jd by the bridging 
element. The smaller Cu-O-Cu angle in HD of 104.98° 



FIG. 5. (a) Predicted evolution of Gu-O-Gu angles as a 

function of pressure for kapellasite, haydeeite, and the hypo¬ 
thetical Cd-kapellasite. (b) Pressure induced changes in the 
two dominant exchange coupling parameters Ji and Jd for 
the three compounds. 
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vs 105.84° in KL results in a larger value for Ji, while 
the additional hopping path via semicore Zn 3d states 
provides a larger J^. The principal discrepancy with 
the experimental values appears to be the overestima¬ 
tion of Jd in HD. Regardless of whether this is an exper¬ 
imental problem {e.g., imperfect samples) or theoretical 
{e.g., overestimation of Mg-0 spa hopping), the trends 
in the dependence of the exchange parameters with re¬ 
spect to both structural and chemical changes are well 
reproduced. Having identified the origin of the behavior 
of Ji and Jd, we propose two recipes for sampling the 
phase diagram. 

The first option is to apply pressure, keeping the chem¬ 
ical compositions. To investigate this avenue, we have 
calculated the structures of KL and HD at experimental 
volumes and at compressions up to 12% (Fig. 5). These 
compressions correspond to pressures of « 7.9 GPa for 
both KL and HD, which is experimentally accessible. In 
both cases, the Cu-O-Cu angle decreases systematically, 
and Ji increases by up to 140% (KL) and 80% (HD), as 
shown in Fig. 5. Thus, by applying pressure we should 
be able to drive KL into the PM regime, and even ap¬ 
proach the boundary with the FM phase, while applying 
pressure to HD will drive it deeper into the FM regime. 

The second option combines both structural and chem¬ 
ical changes. One can preserve homogeneity by substitut¬ 
ing Zn in KL with Cd. Due to the larger ionic radius, this 
substitution may be difficult to realize and might require 
high-pressure synthesis. Indeed, our calculations show 
that while such a compound would be locally stable, the 


Cu-Cd plane would be considerably expanded and the 
Cu-OH-Cu angle would flatten to 112°-113°, which ren¬ 
ders Ji antiferromagnetic [37]. This compound, however, 
would also be highly compressible and would return to a 
structure similar to KL at a pressure of about 20 GPa. 
One expects that Ji should already be ferromagnetic at 
this pressure, while at the same time, we note that since 
the 5d level in Gd lies considerably lower compared to 
the 4:d level in Zn, the Jd is expected to reduce. Both 
conjectures are confirmed by our calculations. Indeed, at 
P = 20 GPa the representative point in the phase dia¬ 
gram appears close to HD at P = 0, while at P = 13.6 
GPa it is close to KL at P = 0 (Fig. 2). Thus, synthesiz¬ 
ing GdGu 3 (OH) 6 Gl 2 and applying external pressure pro¬ 
vides us a vehicle to traverse a vast extent of the phase 
diagram encompassing a large range of Jd, especially deep 
into the paramagnetic phase. 
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I. DENSITY FUNCTIONAL THEORY CALCULATIONS 
A. Prediction of high pressure structures 

We predict high pressure structures using the Vienna ab initio simulation package (vASP) implementation [1, 2] of 
DFT within the generalized gradient approximation (GGA) [3]. We base our calculations on the structures of Krause 
et al. [4] for kapellasite, and of Malcherek et al. [5] for haydeeite. We obtain a P = 0 GPa relaxed structure by fixing 
the experimental volume and relaxing all lattice and internal structural parameters. The DFT calculated pressure at 
this volume, due to the well known underbinding of the GGA exchange correlation functional, is Poffset = 1-54 GPa 
for kapellasite and Poffset = 0.93 GPa for haydeeite. We then proceed to reduce the volume in steps of 3%, relaxing 
all lattice and internal structural parameters. The physical pressure values are obtained by subtracting Poffset from 
the calculated pressure at each reduced volume. 


B. Determination of Hamiltonian parameters 

The calculations were performed with the full potential local orbital (FPLO) basis set [ 6 ] and the GGA functional [3]. 
Total energies for sets of different spin configurations were calculated with GGA+U using the atomic limit double 
counting correction [7], and we employed orthogonal projection of the Cu 3d densities. We would like to note that the 
results are very sensitive to (i) the choice of exchange correlation functional (LDA or GGA), (ii) the double counting 
correction in the LDA+U approach, and (iii) the choice of projection scheme for the Cu 3d orbitals in the FPLO basis. 
In Ref. [ 8 ] the calculations were performed with the LDA exchange correlation functional and the around mean field 
double counting correction was chosen for the LDA+U functional. With this setup the ferromagnetic contribution in 
Ji for kapellasite is strongly suppressed, however, with the same setup the results of our two groups agree [9]. 

We fix Jh = 1 eV and vary [/. Tables I and II show the result of fitting to 9 spin configurations of an orthorhombic 
-v/2 X \/2 X 1 supercell with PI space group, containing 6 inequivalent copper ions. A 6 x 6 x 6 fc mesh was used 
after confirming that an 8 x 8 x 8 A: mesh yields identical results. The error given in Tables I and II reflects the 
statistical error of fitting 9 energies to 5 unknowns (reference energy and four exchange couplings). We have also 
performed calculations with 2x2x1 supercells of kapellasite and haydeeite yielding 46 distinct total energies from 
12 inequivalent copper ions. The two approaches proved to be precisely equivalent so that the computationally less 
demanding smaller supercell was chosen for the final computations. Very small fitting errors indicate that the Cu^+ 
spins are well localized and their interaction is very well described by the chosen Heisenberg model. 

The two significant exchange couplings, Ji and Jd, obtained in this way, are shown visually in Fig. 1 for kapellasite 
and in Fig. 2 for haydeeite. Evolution of couplings with pressure, and with onsite interaction U is very smooth. The 
ferromagnetic nearest neighbour coupling in the kagome plane is larger in haydeeite compared to kapellasite while 
the antiferromagnetic exchange Jg coupling Cu^"*" ions across the hexagons of the kagome lattice is slightly larger in 
kapellasite than in haydeeite. Ji varies more strongly with pressure than Jy, reflecting its sensitivity to the bond 
angle. Increasing the onsite interaction U from U = 6 eV to 1/ = 7 eV and U = 8 eV leads to a smooth decrease in 
absolute value of both exchange couplings. The other two exchange couplings are not shown in the figure; J 2 is very 
small at ~ 0.5 K and insensitive to pressure, while J 4 is nearly exactly zero within computational error. 
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TABLE I: Kapellasite exchange couplings. All structures except the one marked “exp.” are predicted as described In the 
text. 


P (GPa) 

Ji 

J2 

Ji 

Jd 

GGA+U, P = 8 eV, Jn = 1 eV 

0 (exp.) 

-12.59(3) 

-0.55(3) 

-0.17(3) 

16.16(4) 

0 


-0.55(2) 

-0.08(3) 

16.40(4) 

1.6 

-20.78(2) 

-0.50(2) 

-0.04(3) 

17.57(4) 

3.4 

-26.41(2) 

-0.48(2) 

0.01(3) 

18.87(4) 

5.5 

-32.02(2) 

-0.49(2) 

0.04(3) 


7.9 

-37.48(2) 

-0.53(2) 

0.07(2) 

21.80(3) 

GGA+U, P = 7 eV, Jn = 1 eV 

0 

-16.65(4) 

-0.56(4) 

-0.02(4) 

19.89(6) 

1.6 

-22.94(4) 

-0.50(4) 

0.03(4) 

21.33(6) 

3.4 

-29.37(3) 

-0.48(3) 

0.09(4) 

22.95(5) 

5.5 

-35.79(3) 

-0.48(3) 

0.13(4) 

24.71(5) 

7.9 

-42.06(3) 

-0.52(3) 

0.16(4) 

26.58(5) 

GGA+U, P = 6 eV, Jn = 1 eV 

0 

-18.05(6) 

-0.53(6) 

0.10(7) 

24.00(9) 

1.6 

-25.22(6) 

-0.46(6) 

0.17(6) 

25.80(9) 

3.4 

-32.57(5) 

-0.42(5) 

0.23(6) 

27.81(9) 

5.5 

-39.93(5) 

-0.43(5) 

0.28(6) 

29.98(8) 

7.9 

-47.13(6) 

-0.46(6) 

0.28(7) 

32.33(10) 


TABLE II: Haydeeite exchange couplings. All structures except the one marked “exp.” are predicted as described In the text. 


P (GPa) 

Ji 

J2 

Ji 

Jd 

GGA+U, P = 8 eV, J/f = 1 eV 

0 (exp.) 

-21.27(2) 

-0.57(2) 

-0.13(3) 

12.72(4) 

0 

-28.39(2) 

-0.55(2) 

-0.02(3) 

14.06(4) 

1.8 

-33.94(2) 

-0.51(2) 

0.03(2) 

15.42(3) 

3.5 

-39.50(2) 

-0.50(2) 

0.08(2) 

16.97(3) 

5.6 

-45.07(2) 

-0.51(2) 

0.12(2) 

18.71(3) 

7.9 

-50.50(3) 

-0.52(3) 

0.08(3) 

20.65(4) 

GGA+U, P = 7 eV, Jn = 1 eV 

0 

-31.23(3) 

-0.57(3) 

0.05(4) 

16.60(5) 

1.8 

-37.57(3) 

-0.53(3) 

0.11(3) 

18.26(5) 

3.5 

-43.92(3) 

-0.52(3) 

0.16(3) 

20.14(4) 

5.6 

-50.30(3) 

-0.53(3) 

0.21(3) 

22.25(4) 

7.9 

-56.56(2) 

-0.57(2) 

0.25(3) 

24.59(4) 

GGA+U, P = 6 eV, J/f = 1 eV 

0 

-34.15(5) 

-0.57(5) 

0.17(5) 

19.52(7) 

1.8 

-41.38(6) 

-0.52(6) 

0.25(7) 

21.52(9) 

3.5 

-48.64(4) 

-0.50(4) 

0.31(5) 

23.81(6) 

5.6 

-55.94(4) 

-0.52(4) 

0.36(4) 

26.36(6) 

7.9 

-63.13(4) 

-0.57(4) 

0.40(4) 

29.18(6) 


II. ORIGIN OF EXCHANGE PARAMETERS 

Cu2+ is usually a textbook case for the Hubbard model. Cu electrons are strongly localized with hopping t U, 
and the standard superexchange pertnrbation theory works very well. Most experimental papers interpret haydeeite 
(HD) and kapellasite (KL) in terms of this theory, specifically, in terms of the Goodenough-Kanamori-Anderson 
rules, which stipulate that at the bond angle (j) = 90°, the superexchange process Cu-O-Cu is completely suppressed 
and the only remaining interaction is ferromagnetic and generated by the Hund’s rule coupling on oxygen (the fact 
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TABLE III: Hypothetical Cd-kapellasite exchange couplings. All structures are predicted as described In the text. 


P (GPa) 

Ji 

J 2 

J 4 

Jd 

GGA+U, G = 8 eV, J/f = 1 eV 

0.6 

36.96(5) 

-0.10(3) 

-0.46(5) 

7.15(7) 

3.2 

23.82(4) 

-0.09(2) 

-0.33(5) 

8.04(6) 

4.8 

17.13(4) 

-0.08(2) 

-0.27(4) 

8.55(6) 

7.6 

6.98(3) 

-0.08(2) 

-0.20(4) 

9.38(5) 

9.8 

0.25(3) 

-0.08(2) 

-0.15(4) 

9.99(5) 

13.6 

-9.95(3) 

-0.08(2) 

-0.11(3) 

11.00(5) 

20 

-23.31(3) 

-0.09(2) 

-0.11(3) 

12.54(4) 

GGA+U, U = 7 eV, Jh = ley 

0.6 

42.37(7) 

-0.10(4) 

-0.46(8) 

8.77(11) 

3.2 

27.39(6) 

-0.08(3) 

-0.31(7) 

9.90(9) 

4.8 

19.74(5) 

-0.08(3) 

-0.25(6) 

10.53(9) 

7.6 

8.10(5) 

-0.07(3) 

-0.15(6) 

11.58(8) 

9.8 

0.37(5) 

-0.07(3) 

-0.11(5) 

12.34(7) 

13.6 

-11.39(4) 

-0.08(2) 

-0.06(5) 

13.62(7) 

20 

-26.83(4) 

-0.09(2) 

-0.08(4) 

15.53(6) 

GGA+U, G = 7 eV, Jff = 1 eV 

0.6 

48.46(10) 

-0.09(5) 

-0.41(11) 

10.69(16) 

3.2 

31.40(9) 

-0.074(5) 

-0.23(10) 

12.13(14) 

4.8 

22.67(7) 

-0.07(4) 

-0.13(8) 

12.93(11) 

9.8 

0.42(7) 

-0.06(4) 

0.00(8) 

15.20(11) 

13.6 

-13.14(7) 

-0.07(4) 

0.05(8) 

16.78(11) 

20 

-31.05(6) 

-0.08(3) 

-0.01(6) 

19.14(10) 


that in these materials the angle is actually rather far from 90° is usually swept under the rug). This coupling is 
proportional to (where tpd is the hopping and Ap^ is the charge transfer energy). Note that the usual 

antiferromagnetic superexchange is proportional to which is larger by a factor A^JU Jq for A'^JU Jq > 1- 

Thus, near-cancellation of the assisted hopping via different O orbitals is a necessary requirement for ferromagnetic 
superexchange. There are two ways to look at this situation. One can assume a coordinate system where x and y 
are along the Cu-0 and 0-Cu bonds in the trimer, in which case the “left” Cu electron can only hop to the O 
orbital, and the “right” one only to the O py one. Since these are mutually orthogonal, there is no net assisted Cu-Cu 
hopping. Another way is to assume the system where x is along the Cu-Cu bond and y is perpendicular to it. In 
this case, there is assisted hopping via each of the orbitals, but the two contributions are equal in magnitude and 
opposite in sign, so they cancel out. As the angle becomes larger than 90°, the Cu-pa;-Cu hopping becomes larger 
than Cu-pp-Cu one and cancellation becomes incomplete. In typical transition metal oxides if the angle is larger than 
« 95° this effect is significant. 

However, our calculations as well as the experimental results, indicate that the crossover from the antiferromagnetic 
to ferromagnetic Ji in these materials occurs at the bond angle of 108°—109°. Moreover, as discussed in Ref. [8] and 
confirmed in our calculations, shifting H away from O makes the interaction at a given angle less ferromagnetic, 
suggesting an important role of H orbitals. 

In the literature one can find two seemingly contradicting predictions regarding adding to the standard 4-orbitals 
(two Cu d and two O d) another bridging orbital. On one hand, Geertsma and Khomskii considered a (j) = 90° case 
[10] when O p orbitals are bound to another orbital, such as in our case hydrogen, in their case germanium, or, in 
the simplest case with the same symmetry, to O s. In that case the py orbital couples with H, but Px does not. As 
a result, the cancellation described above is incomplete and there is some residual antiferromagnetic superexchange, 
which can either overcome the weak ferromagnetic coupling, as in Ref. [10], or considerably reduce it. Furthermore, 
this theory predicts that at </) = 90° moving H toward O must enhance the tendency to antiferromagnetism. 

On the other hand. Hay et al. ]11] considered a similar problem (using O s and the additional bridging orbital, 
which is not important) as a function of (j) and concluded that the full cancellation between all Cu-O-Cu hopping 
processes occurs not at 0 = 90°, but at a larger [for the Cu 2 (OH) 2 Cl 4 ~ radical they found numerically (j)c « 108°, 
but the actual value must be system-dependent, and, in particularly, may be smaller for our materials]. This is not 
in contradiction with Ref. ]10]; indeed Hay et al. explicitly emphasize that it is only at 0c that the antiferromagnetic 
exchange completely cancels out; for both larger and smaller angles there is residual antiferromagnetic coupling, in 
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kapellasite 



FIG. 1: Exchange couplings of kapellasite. The two significant Heisenberg Hamiltonian parameters for kapellasite are given 
as function of applied pressure and for three different values of the onsite interaction strength U. 


accord with the Geertsma and Khomskii theory. However, at </> = 90° the Cu-py-Cu hopping is the stronger the 
shorter the 0-H bond is, and to restore the balance one needs to enhance Cu-p^^-Cu hopping by flattening cj). So, at 
least qualitatively (quantitative analysis is underway and will be published elsewhere). Hay et al. theory seems to 
explain why the large bond angles in HD and KL nontheless generate net ferromagnetic coupling. 

Yet another paradox manifests itself in comparing the Jd values for HD and KL. It is generally believed that this 
exchange is provided by hopping from Cu {via O) to Zn or Mg. However, as discussed below, this may be not the 
case. Indeed, Cu-(Zn,Mg) hopping can proceed either to the s orbitals of Zn or Mg, or to d orbitals of Zn. The 5s 
orbitals of Zn are considerably more diffuse than the 3s orbitals of Mg, and much closer to the Fermi level. Since Jd 
in the simple Hubbard theory should be proportional to the square of the overlap integrals and inversely proportional 
to the energy separation, we estimate that Jd in HD should be nearly an order of magnitude smaller than in KL, 
even without taking into account an additional hopping channel via 3d orbitals of Zn. In the experiment, as well as 
in the calculations, the difference is less than a factor of two. A closer look, however, shows that this naive reasoning 
doesn’t apply here. Indeed, an inspection of the crystal structure shows that the Cu-02-(Zn,Mg)-02-Cu multimer is 
planar, as shown in Fig. 4. The bond angle is about 100° considerably closer to 90° that the Cu-O-Cu angle. Had it 
been exactly 90°, the hopping from the Cu orbital, where the actual hole resides, to the (Zn,Mg) s orbital 

would be zero by parity (see Fig. 4(b)). For 100°, it should be strongly suppressed. For Zn, there is an additional 
option of hopping via 3d orbitals. However, as Fig. 4(c) shows, dpa — pda hopping is again forbidden by symmetry. 
The dpa — pdir hopping is allowed, but the two oxygens are antiphase and cancel each other (see Fig. 4(d)). 

Keeping in mind the H effect, one may conjecture that 100° may be rather close to full cancellation. Our numerical 
results provide for the effective (after integrating out O) Cu(d)—Zn(d) hopping a value of about r fa 40 meV, which 
is nearly an order of magnitude smaller that for the Cu-Cu one. Thus, we conclude that the mechanism providing 
for sizeable Jd in HD and KL, as compared to herbertsmithite, is more complex than it has been anticipated so 
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haydeeite 



FIG. 2: Exchange couplings of haydeeite. The two significant Heisenberg Hamiltonian parameters for haydeeite are given 
as function of applied pressure and for three different values of the onsite interaction strength U. 


far, and does not reduce to just providing a hopping path via a central atom. This makes it very difficult, if at 
all possible, to anticipate the trends in the exchange parameters after material modifications and makes ab initio 
calculations indispensable. As a last remark, it should be mentioned that for a systematic investigation of the effects 
of, for instance, substitution of the central atom, experimental estimates are not necessarily superior to the theoretical 
estimates, due to the problems related to site disorder in the samples. 
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FIG. 3: Exchange couplings of hypothetical Cd-kapellasite. The two significant Heisenberg Hamiltonian parameters for 
Cd-kapellasite are given as function of applied pressure and for three different values of the onsite interaction strength U. 
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FIG. 4: Orbitals and hopping pathways. 


